﻿Imports System.Math


<Serializable()> Public Class tunnelArc
    Public Radius As Double
    Public Center As New TSPCoor(0, 0)
    Public StartAngle As Double
    Public EndAngle As Double



    '测量点坐标(Xi,Yi),i=1,2,3,...,N
    '以测量点到最小二乘半径的平方的残差的平方和为最小
    '圆心(Xc,Yc)，半径Rc
    '圆拟合
    Function CircleDefineByPoints(ByRef Pnts As List(Of SectionPoint)) As Boolean
        Dim mPnts As New List(Of SectionPoint)
        For Each pt As SectionPoint In Pnts
            If pt.IsValid Then
                mPnts.Add(pt)
            End If
        Next

        Dim PntCnt As Integer
        PntCnt = mPnts.Count
        If PntCnt < 3 Then Exit Function

        Dim N As Integer : N = PntCnt
        Dim S As Double : S = Tot_One(mPnts, TSPCoor.CoorVH.Vertical)
        Dim T As Double : T = Tot_One(mPnts, TSPCoor.CoorVH.Horizontal)
        Dim H As Double : H = Tot_Two(mPnts, TSPCoor.CoorVH.Vertical, TSPCoor.CoorVH.Horizontal)
        Dim P As Double : P = Tot_Two(mPnts, TSPCoor.CoorVH.Vertical, TSPCoor.CoorVH.Vertical)
        Dim Q As Double : Q = Tot_Two(mPnts, TSPCoor.CoorVH.Horizontal, TSPCoor.CoorVH.Horizontal)
        Dim U As Double : U = Tot_Three(mPnts, TSPCoor.CoorVH.Vertical, TSPCoor.CoorVH.Vertical, TSPCoor.CoorVH.Horizontal)
        Dim V As Double : V = Tot_Three(mPnts, TSPCoor.CoorVH.Vertical, TSPCoor.CoorVH.Horizontal, TSPCoor.CoorVH.Horizontal)
        Dim W As Double : W = Tot_Three(mPnts, TSPCoor.CoorVH.Vertical, TSPCoor.CoorVH.Vertical, TSPCoor.CoorVH.Vertical)
        Dim Z As Double : Z = Tot_Three(mPnts, TSPCoor.CoorVH.Horizontal, TSPCoor.CoorVH.Horizontal, TSPCoor.CoorVH.Horizontal)

        Dim TT As Double
        Dim Xc As Double, Yc As Double, Rc As Double
        TT = (P + P - 2 * S * S / N) * (Q - T * T / N) - (H + H - 2 * S * T / N) * (H - S * T / N)
        If TT = 0 Then
            MsgBox("Err:圆拟合时除数为0")
            CircleDefineByPoints = False
        Else
            Xc = (Q - T * T / N) * (W + V - S * (P + Q) / N) - (H - S * T / N) * (Z + U - T * (P + Q) / N)
            Xc = Xc / TT
            Yc = (P - S * S / N) * (Z + U - T * (P + Q) / N) - (H - S * T / N) * (W + V - S * (P + Q) / N)
            Yc = Yc / TT
            Rc = Math.Sqrt(Xc * Xc + Yc * Yc - 2 * S * Xc / N - 2 * T * Yc / N + (P + Q) / N)

            Center.Coor_V = Xc
            Center.Coor_H = Yc
            Radius = Rc
            CircleDefineByPoints = True

        End If
    End Function

    '求各点与拟合圆的偏差
    Public Sub Cal_Err_From_Cir(ByRef Sections As TunnelSectionPoints, ByRef Errs() As Double, ByRef StandErr As Double)
        If Sections.Points.Count = 0 Then
            Return
        End If
        Dim Pnts As List(Of SectionPoint) = Sections.Points
        ReDim Errs(Pnts.Count - 1)
        Dim ValidPointCount As Single = 1.0E-30
        Dim mErr As Double = 0

        Dim I As Integer
        For I = 0 To Pnts.Count - 1
            Dim vP As SectionPoint = Pnts(I)
            Dim dv As Double = vP.Coor.Coor_V - Center.Coor_V
            Dim dh As Double = vP.Coor.Coor_H - Center.Coor_H

            Errs(I) = Math.Sqrt(dv * dv + dh * dh) - Radius
            Sections.Points(I).Err_Simulation = Errs(I)
            If vP.IsValid Then
                ValidPointCount += 1
                mErr += Errs(I) * Errs(I)
            End If
        Next
        StandErr = Math.Abs(mErr / (Pnts.Count))
        Sections.StdErr = StandErr
        If (Radius > 0) Then

            Sections.center = Center

            Sections.Radius = Radius

        End If
        Sections.start_hv = Sections.Points(0).AzimFrom(Sections.center)
        Sections.end_hv = Sections.Points(Sections.Points.Count - 1).AzimFrom(Sections.center)
        Sections.ArdLength = Radius * (Sections.end_hv - Sections.start_hv)



    End Sub





    '计算三点确定的圆形的圆心和半径
    Function ArcDefBy3P(ByVal pt1 As TSPCoor, ByVal pt2 As TSPCoor, ByVal pt3 As TSPCoor) As Boolean
        Dim xysm, xyse, xy As Double
        xy = pt1.Coor_V ^ 2 + pt1.Coor_H ^ 2
        xyse = xy - pt3.Coor_V ^ 2 - pt3.Coor_H ^ 2
        xysm = xy - pt2.Coor_V ^ 2 - pt2.Coor_H ^ 2
        xy = (pt1.Coor_V - pt2.Coor_V) * (pt1.Coor_H - pt3.Coor_H) - (pt1.Coor_V - pt3.Coor_V) * (pt1.Coor_H - pt2.Coor_H)
        '判断参数有效性
        If Math.Abs(xy) < 0.000001 Then
            'MsgBox("所输入的参数无法创建圆形！")
            Return False
            Exit Function
        End If
        '获得圆心和半径
        Center.Coor_V = (xysm * (pt1.Coor_H - pt3.Coor_H) - xyse * (pt1.Coor_H - pt2.Coor_H)) / (2 * xy)
        Center.Coor_H = (xyse * (pt1.Coor_V - pt2.Coor_V) - xysm * (pt1.Coor_V - pt3.Coor_V)) / (2 * xy)
        Radius = Sqrt((pt1.Coor_V - Center.Coor_V) ^ 2 + (pt1.Coor_H - Center.Coor_H) ^ 2)
        Radius = Sqrt((pt2.Coor_V - Center.Coor_V) ^ 2 + (pt2.Coor_H - Center.Coor_H) ^ 2)
        Radius = Sqrt((pt3.Coor_V - Center.Coor_V) ^ 2 + (pt3.Coor_H - Center.Coor_H) ^ 2)
        Return True
    End Function

    'Dim xysm, xyse, xy As Double

    'xy = pt1(0) ^ 2 + pt1(1) ^ 2
    'xyse = xy - pt3(0) ^ 2 - pt3(1) ^ 2
    'xysm = xy - pt2(0) ^ 2 - pt2(1) ^ 2
    'xy = (pt1(0) - pt2(0)) * (pt1(1) - pt3(1)) - (pt1(0) - pt3(0)) * (pt1(1) - pt2(1))
    ''判断参数有效性
    'If Abs(xy) < 0.000001 Then
    '    MsgBox "所输入的参数无法创建圆形！"
    '    GetCenOf3Pt = False
    '    Exit Function
    'End If
    ''获得圆心和半径
    'PtCen(0) = (xysm * (pt1(1) - pt3(1)) - xyse * (pt1(1) - pt2(1))) / (2 * xy)
    'PtCen(1) = (xyse * (pt1(0) - pt2(0)) - xysm * (pt1(0) - pt3(0))) / (2 * xy)
    'PtCen(2) = 0
    'Radius = Sqr((pt1(0) - PtCen(0)) * (pt1(0) - PtCen(0)) + (pt1(1) - PtCen(1)) * (pt1(1) - PtCen(1)))
    'GetCenOf3Pt = True




    '计算二点及半径确定的圆弧圆心
    Function ArcDef_By_2P_And_R(ByVal P1 As TSPCoor, ByVal P2 As TSPCoor, ByVal Radius As Double, ByVal CenUp As Boolean) As Boolean
        If P1 Is Nothing Then
            Throw New ApplicationException("Exception Occured")
        End If
        Dim sVct As New SvVector(P1, P2)
        If sVct.Distance < Radius * 2 Then
            Dim pMid As New TSPCoor(P1.Coor_V / 2 + P2.Coor_V / 2, P1.Coor_H / 2 + P2.Coor_H / 2)
            Dim sDst As Double = Math.Sqrt(Radius * Radius - sVct.Distance * sVct.Distance / 4)

            Dim sAzim As Double = sVct.Azimuth - Math.PI / 2
            If CenUp Then
                If sDst * Math.Cos(sAzim) < 0 Then sAzim += Math.PI
            Else
                If sDst * Math.Cos(sAzim) > 0 Then sAzim += Math.PI
            End If

            Center.Coor_V = pMid.Coor_V + sDst * Math.Cos(sAzim)
            Center.Coor_H = pMid.Coor_H + sDst * Math.Sin(sAzim)

            Me.Radius = Radius
        Else
            Return False
        End If

        Return True
    End Function



    '三次方求和
    Private Function Tot_Three(ByRef mPnts As List(Of SectionPoint), ByVal XorY1 As TSPCoor.CoorVH, ByVal XorY2 As TSPCoor.CoorVH, ByVal XorY3 As TSPCoor.CoorVH)
        Dim Tot As Double
        For Each p As SectionPoint In mPnts

            Dim tmp1 As Double = IIf(XorY1 = TSPCoor.CoorVH.Vertical, p.Coor.Coor_V, p.Coor.Coor_H)
            Dim tmp2 As Double = IIf(XorY2 = TSPCoor.CoorVH.Vertical, p.Coor.Coor_V, p.Coor.Coor_H)
            Dim tmp3 As Double = IIf(XorY3 = TSPCoor.CoorVH.Vertical, p.Coor.Coor_V, p.Coor.Coor_H)
            Tot = Tot + tmp1 * tmp2 * tmp3

        Next
        Tot_Three = Tot
    End Function

    'X*Y /X*X /Y*Y 求和
    Private Function Tot_Two(ByRef mPnts As List(Of SectionPoint), ByVal XorY1 As TSPCoor.CoorVH, ByVal XorY2 As TSPCoor.CoorVH)

        Dim Tot As Double
        For Each p As SectionPoint In mPnts

            Dim tmp1 As Double = IIf(XorY1 = TSPCoor.CoorVH.Vertical, p.Coor.Coor_V, p.Coor.Coor_H)
            Dim tmp2 As Double = IIf(XorY2 = TSPCoor.CoorVH.Vertical, p.Coor.Coor_V, p.Coor.Coor_H)
            Tot = Tot + tmp1 * tmp2
        Next
        Tot_Two = Tot
    End Function

    'X 或 Y 求和
    Private Function Tot_One(ByRef mPnts As List(Of SectionPoint), ByVal XorY As TSPCoor.CoorVH) As Double
        Dim Tot As Double

        For Each p As SectionPoint In mPnts
            Dim tmp1 As Double = IIf(XorY = TSPCoor.CoorVH.Vertical, p.Coor.Coor_V, p.Coor.Coor_H)
            Tot = Tot + tmp1
        Next
        Tot_One = Tot
    End Function

End Class
